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Multiscale simulation methods have been developed based on the local stress sampling strategy 
and applied to three flow problems with different difficulty levels: (a) general flow problems of simple 
fluids, (b) parallel (one-dimensional) flow problems of polymeric liquids, and (c) general (two- or 
three-dimensional) flow problems of polymeric liquids. In our multiscale methods, the local stress 
of each fluid element is calculated directly by performing microscopic or mesoscopic simulations 
according to the local flow quantities instead of using any constitutive relations. For simple fluids 
(a), such as the Lenard- Jones liquid, a multiscale method combining MD and CFD simulations is 
developed based on the local equilibrium assumption without memories of the flow history. The 
results of the multiscale simulations are compared with the corresponding results of CFD with 
or without thermal fluctuations. The detailed properties of fluctuations arising in the multiscale 
simulations are also investigated. For polymeric liquids in parallel flows (b), the multiscale method 
is extended to take into account the memory effects that arise in hydrodynamic stress due to the slow 
relaxation of polymer-chain conformations. The memory of polymer dynamics on each fluid element 
is thus resolved by performing MD simulations in which cells are fixed at the mesh nodes of the CFD 
simulations. The complicated viscoelastic flow behaviours of a polymeric liquid confined between 
oscillating plates are simulated using the multiscale method. For general (two- or three-dimensional) 
flow problems of polymeric liquids (c), it is necessary to trace the history of microscopic information 
such as polymer-chain conformation, which carries the memories of past flow history, along the 
streamline of each fluid element. A Lagrangian-based CFD is thus implemented to correctly advect 
the polymer-chain conformation consistently with the flow. On each fluid element, coarse-grained 
polymer simulations are carried out to consider the dynamics of entangled polymer chains that show 
extremely slow relaxation compared to microscopic time scales. This method is successfully applied 
to simulate a flow passing through a cylindrical obstacle. 

PACS numbers: 31.15.xv 46.15.-x 

Keywords: multiscale, simulation, modeling, polymeric liquids, complex fluids, softmatters, viscoelastic flu- 
ids, memory effect 
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I. INTRODUCTION 



The prediction of flow behaviours of polymeric liquids 
using computer simulations is a challenging theme in var- 
ious fields of science and engineering including physical 
science, materials science, and mechanical and chemical 
engineering. Polymeric liquids are known to exhibit pe- 
culiar flow behaviours that are related to the microscale 
dynamics of their polymer chains; these dynamics affect 
the viscoelasticity, shear thinning/thickening behaviours, 
and flow-induced phase transitions of these liquids^ The 
characteristic times of the microscale dynamics of poly- 
mer chains tend to be very long and are often compara- 
ble to the time scales of macroscale fluid motions. Thus, 
for these compounds, one cannot separate the microscale 
and macroscale dynamics in the temporal domain, even 
though a large gap exists between those length and time 
scales; one also cannot make the assumption usually 
made for simple fluids under flow that the fluid elements 
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are in local equilibrium at each instance. Therefore, one 
must trace the "memory" or "history" of the deforma- 
tions of each fluid element along its streamline. This 
coupling between microscale and macroscale dynamics 
hinders the simulation of polymeric liquids. 

In this paper, we present a detailed explanation of a 
multiscale method that bridges the hydrodynamic mo- 
tions of fluids using computational fluid dynamics (CFD) 
and the microscopic [or mesoscopic] dynamics of poly- 
mer configurations using molecular-dynamics (MD) [or 
coarse-grained (CG)2~— ] simulations. The concept of 
bridging microscale and macroscale dynamics is also im- 
portant for other flow problems of softmatters with com- 
plex internal degrees of freedom (e.g., colloidal disper- 
sion, liquid crystal, and glass). The basic idea of our 
multiscale method can be applied to those softmatter 
flows. 

In non-Newtonian fluid dynamics, many novel CFD 
schemes have been proposed for the viscoelastic flows of 
polymeric liquids»2r— In CFD methods, a model consti- 
tutive equation is used to determine the local stresses 
at each instant from the history of previous velocity 
fields. For dense polymeric liquids (melts), however, the 
detailed form of the constitutive relations is so compli- 
cated that they are generally unknown. 1 - Thus, the usual 
CFD methods, which require constitutive equations, are 
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FIG. 1: Schematic illustrations of our multiscale methods bridging CFD and MD [or CG] simulations. The developments 
were carried out in a step-by-step manner starting from the bridging method for simple fluids (a) and proceeding to a method 
applicable to general flow problems of complex fluids. The small cubic boxes placed in the fluid systems represent MD or CG 
simulation cells within which the local stress of each fluid element is sampled. Case (a) represents a cavity flow of a simple 
fluid for which the local stress a of the fluid at a time t and a position r is given by a function of the local deformation rate 
7 at the same time and position, a(r,t) — /(7(r,£)). Case (b) represents a polymeric liquid subject to an oscillatory shear 
flow where the flow is parallel to the a;-direction and the velocity gradient exists only in the y-direction. Here, the local stress 
is a function of the history of the velocity gradient in the past t' < t at the position y, a(y,t) = F[-y(y,t'); t' < t]. Case (c) 
represents a polymeric liquid passing through a cylindrical obstacle. Here, the local stress is a function of the history of the 
velocity gradient in the past t' < t along the streamline R(t') of the fluid element, a(r,t) — F[j(R(t'),t'); t' <t]. The details 
of our multiscale methods for cases (a), (b), and (c) are given in Sees. II, III, and IV, respectively. 



usually not straightforwardly applicable to the compli- 
cated flow problems of polymeric liquids. Instead, mi- 
croscopic simulations such as MD and CG simulations 
are often used to investigate the rheological properties 
of such materials. The microscopic simulations are usu- 
ally performed only for a tiny piece of the material in 
the equilibrium or non-equilibrium state under uniform 
external fields of shear velocity, temperature gradient, 
and electric field. Although the microscopic simulations 
are even applicable to macroscopic flow behaviours, the 
drawback of this type of simulation is the enormous com- 
putation time required. Thus, for problems that concern 
large-scale and long-time fluid motions far beyond the 
molecular scale, which are commonly encountered in en- 
gineering, fully microscopic simulations are difficult from 
a practical standpoint. To overcome the weaknesses of 
the individual methods, we have developed a multiscale 
simulation method that combines CFD and MD [or CG] 
simulations fibr- 
in our multiscale methods, macroscopic flow behaviour 
is calculated by a CFD solver; however, instead of us- 
ing any constitutive equations, a local stress is calcu- 
lated using the MD [or CG] simulation associated with 
a local fluid element according to the local flow variable 
obtained by the CFD calculation. The multiscale meth- 
ods are applied to three flow problems with different lev- 
els of difficulty: (a) general flow problems of simple flu- 



ids, (b) parallel (one-dimensional) flow problems of poly- 
meric liquids, and (c) general (two- or three-dimensional) 
flow problems of polymeric liquids, as schematically illus- 
trated in Fig. [T] 

For the simple fluids shown in Fig. [T] (a), the local 
stress a of the fluid at a time t and a position r is given 
by a function of the local deformation rate j at the same 
time and position, 

a(r,t) = f(j(r,t)). (1) 

We proposed a multiscale method that combines CFD 
and MD based on the local equilibrium assumption, in 
which a local stress immediately attains a steady state 
after a short transient time during which a strain rate is 
subjected to the fluid element. In this method, a lattice- 
mesh-based CFD scheme is used at the macroscopic level, 
and the non-equilibrium MD simulations are performed 
in small MD cells associated with each lattice node of the 
CFD to generate a local stress according to a local strain 
rate. 

The multiscale method combining CFD and MD is ex- 
tended in a straightforward fashion to the polymeric flows 
in the one-dimensional geometry shown in Fig. Q] (b), in 
which the macroscopic quantities, e.g., velocity, temper- 
ature, and stress, are uniform in the flow direction par- 
allel to the plates. This situation allows us to neglect 
the advection of memory on a fluid element along the 



3 



streamline, so the local stress of a fluid element at t and 
y is a functional of the history of the velocity gradient in 
the past t' < t at the same position y, 

a(y,t) = F[j(y,t'); t' < t}. (2) 

For general two- or three-dimensional flow problems of 
polymeric liquids, shown in Fig.[T](c), one must consider 
the advection of polymer-chain conformations, which car- 
ries memory effects on a fluid element along the stream- 
line. The local stress is thus a functional of the history of 
the velocity gradient in the past t' <t along the stream- 
line R(t') of the fluid element, 

a(r,t)=F[j(R(t'),t'); t'<t]. (3) 

To meet this requirement, Lagrangian fluid dynamics are 
implemented for the CFD calculation to trace the advec- 
tion of a fluid element that contains the memory of the 
configuration of the polymer chains. The local stresses 
on each fluid particle are calculated using coarse-grained 
polymer simulations in which the dynamics of entangled 
polymer chains are calculated. 

The idea of using multiscale modelling to calculate 
the local stress for the fluid solver using the micro- 
scopic simulation instead of any constitutive relations 
was first proposed for polymeric liquids by Laso and 
Ottingepi^SS, who presented the CONNFFESSIT ap- 
proach. The multiscale method was also proposed by 
E and Enquist^i who presented the heterogeneous mul- 
tiscale method (HMM) as a general methodology for the 
efficient numerical computation of problems with multi- 
scale characteristics. HMM has been applied to the sim- 
ulation of complex fluid o 22 ' 23 but completely neglects the 
advection of memory. The equation-free multiscale com- 
putation proposed by Kevrekidis et al. is based on a sim- 
ilar idea and has been applied to various problemsi 24 i 25 
The basic idea of our multiscale modelling is the same 
as those earlier proposed. This type of bridging method 
has also been developed recently by several researchers. 
De et al. have proposed the scale-bridging method, which 
can correctly reproduce the memory effect of a polymeric 
liquid, and demonstrated the non-linear viscoelastic be- 
haviour of a polymeric liquid between oscillating plates. 25 
The methodology of the present multiscale simulations 
for one-dimensional flows of polymeric liquids is the same 
as that used in the scale-bridging method. Kessler et al. 
have recently developed a multiscale simulation for rar- 
efied gas flows based on a similar idea^l 

In the following, the bridging method of MD and CFD 
simulations for simple fluids is presented in Sec. II. In 
this section, the results of the 2-D cavity flows are com- 
pared with those of Newtonian fluids to demonstrate the 
validity of our multiscale simulation method. A special 
focus is placed on the fluctuation arising in our multiscale 
method. We carry out spectral analysis of the fluctua- 
tions and compare the results with those obtained us- 
ing fluctuating hydrodynamics. In Sec. Ill, the multi- 
scale method of MD and CFD simulations is extended 




FIG. 2: The stress relaxation function G(t) for the LJ liquid 
(dashed line) and model polymer melt (solid line). 

to the one-dimensional flows of polymeric liquids con- 
fined in parallel plates. The flow behaviours of a glassy 
polymer melt in the oscillating plates are calculated, the 
local rheological properties of the polymer melt in the 
rapidly oscillating plates are investigated, and the results 
are compared with the analytical solution for an infinites- 
imally small strain to demonstrate the validity of the 
multiscale method. In Sec. IV, the multiscale method 
of Lagrangian dynamics and CG polymer-dynamics sim- 
ulations is presented for a two-dimensional flow problem 
of polymeric liquids. A transient flow passing a circular 
object is demonstrated. Summaries of each section and 
future perspectives are given in Sec. V. 



II. SIMPLE FLUID 

We consider a simple liquid composed of particles in- 
teracting via the repulsive part of the Lennard- Jones (L J) 
potential, 

_ f 4 £ [(a/r) 12 - [a/rf] + e (r < 2 1 '**), 
U ^-{ (r ' 

(4) 

where a and e are the length and energy units of the LJ 
potential, respectively. In this section, we measure space 
and time in the units of a and tq = ^Jmo 2 Je, respec- 
tively. The temperature T is measured in the unit of 
e/ks- Here, m and kg are the mass of the LJ particle 
and the Boltzmann constant, respectively. The temper- 
ature T and density p of the liquid are assumed to be 
uniform and fixed as T = 1.0 and p — 0.8. 

At this density and temperature, the LJ liquid does 
not have a long-time memory. That is, the shear stress 
depends only on the instantaneous strain rate, regardless 
of the history of the previous strain rates experienced 
by the fluid element. In Fig. [51 we show the stress- 
relaxation function G(t) of the present LJ liquid and of 



4 



a model polymeric liquid (the latter will be discussed in 
the next section). The stress-relaxation function G(t) is 
calculated as 



G(t) = (n xy (t + t )mo))/k B TV 7 



(5) 



where H xy is the space integral of the microscopic stress 
tensor in the volume V. It is seen that the stress relax- 
ation of the LJ liquid rapidly decreases and that the time 
correlation of stress almost disappears in t -C 1. The 
relaxation time tlj of the LJ liquid may be estimated 
as t lj = 0.067 with the definition G(t lj )/G(0) = e _1 . 
The time scale of temporal variations in the macroscopic 
flows may be much larger than the stress-relaxation time 
of the present LJ liquid. Thus, for the macroscopic flow 
behaviours of simple fluids, one can ignore the memory 
effect in a local stress due to the history of previous flow 
velocities. That is, we can assume the local equilibrium 
states at any instant in the macroscopic flows. 

In this section, we present a model for multiscale sim- 
ulations of MD and CFD for simple fluids without mem- 
ory effects. The bridging scheme of MD and CFD is 
based on the local equilibrium assumption of the macro- 
scopic quantities. The multiscale simulations are per- 
formed for the driven cavity flows of the simple LJ liquid, 
and the results are compared with those of the Newto- 
nian fluid. Special attention is given to the efficiency of 
the present multiscale method and to the noise arising in 
this method. 



A. Multiscale Model 

Incompressible flows with uniform density po and tem- 
perature To are described by the following equations: 



dv a 
dx a 



= 0, 



dt 



- v 



dv a 1 dP a f 



dxr- 



Po dxt 



9o 



(6) 



(7) 



where x a is the Cartesian coordinate system, t the time, 
v a the velocity, po the density, P a p the stress tensor, and 
g a the external force per unit mass. Throughout this 
work, the subscripts a, j3, and 7 represent the index in 
Cartesian coordinates, i.e., {a, /?, 7} = {x, y, z}, and the 
summation convention is used. The stress tensor P a p is 
written in the form 



(8) 



where p is the pressure and 8 a p is the Kronecker delta. 
We may assume that T a p is symmetric, T a p — Tp a , and 
traceless, T QQ =0, for isotropic simple fluids^ To solve 
the above equations, one needs a constitutive relation for 
the stress tensor T a p. In our multiscale method, instead 
of using any explicit formulas such as the Newtonian con- 
stitutive relation, T a p is computed directly by MD sim- 
ulations. 



1. CFD Scheme 

We use a lattice-mesh-based finite-volume method 
with a staggered arrangement for vector and scalar 
quantities^ (Fig. [3]). The control volume for a vector 
quantity is a unit square surrounded by dashed lines, and 
that for a scalar quantity is a unit square surrounded by 
solid lines. Equations. © and are discretised by in- 
tegrating the quantities on each control volume. For nu- 
merical time integrations, we use the fourth-order Runge- 
Kutta method, in which a single physical time step At 
is divided into four sub-steps. The time evolution of a 
quantity <fr, which is to be determined by the equation 
d(f>/dt—f(t,(j)), is written as 



At 
At 



f(tn,<t> n ), 



>J n+l - 
in+1 



+Atf(t n+i ,r: +i ), 

At r 



(9a) 

(9b) 
(9c) 
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(9d) 

The time evolution of the fluid velocity v is computed 
by the above set of equations. On the other hand, the 
pressure p is determined so that the fluid velocity satisfies 
the incompressible condition © at each sub-step. The 
procedure at each sub-step is written as 



p = p + ip, 

V = v — tV'0, 

AV> = -Vv, 

T 



(10a) 
(10b) 

(10c) 



where p is the pressure obtained at the previous sub-step, 
v is the velocity obtained by solving equation (|9|) at the 
present sub-step, r is the time increment of the sub-step, 
and ip is the correction term to obtain the divergence- 
free velocities. The remaining three components of the 
stress tensor, T a p, are to be computed directly by MD 
simulations. The details of the method are described in 
the next subsection. Note that the calculation of T a p is 
carried out at each sub-step of Eq. ([9]) . 



2. Local MD sam 



We compute the local stresses by MD simulations ac- 
cording to the local strain rates, rather than the local 
flow velocities themselves. A schematic diagram of the 
method is depicted in Fig. [3J At the CFD level, the local 
strain rate tensor E a p is defined as 



E, 



a/3 



1 / dVa_ 

2 \dx p 



dx a 



(11) 
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FIG. 3: Schematic diagram of the multiscale simulation method for simple fluids, (a) Staggered arrangement of the velocity v 
and stress tensor P a p on a lattice-mesh grid of the CFD calculation. CFD simulations are performed in a reference coordinate 
(x,y,z), while MD simulations are performed in a rotated coordinate (x',y',z') so that the diagonal components of E' a p all 
become zero using the procedure described in Sec II B. The CFD system is discretised into cubic subsystems whose side length 
is Ax. Each subsystem is associated with an MD cell whose side length is Imd, with the Lees-Edward periodic boundary 
condition under shear deformation. Note that three-dimensional MD simulations are used at the microscopic level even for the 
problems for which one- or two-dimensional analysis is applied at the macroscopic level, (b) A schematic time evolution of our 
multi-scale method. The CFD simulation proceeds with a time step of At, and the MD simulation is carried out for a lapse of 
time £md only to sample the local stress T' a p at each node point and time step of CFD. 



where the incompressible condition, E aa =0, is to be sat- 
isfied. We can now define a rotation matrix by which 
the strain-rate tensor E a p is transformed to 



E' = eEe T = 



E 'xy E 'xz 
E 'yx E 'yz 
E zx E 'zy 



(12) 



where the diagonal components all vanish. This transfor- 
mation enables us to perform MD simulations with the 
Lees-Edwards periodic boundary condition. 

We note that the Lees- Edwards periodic boundary con- 
dition cannot create a flow field for an arbitrary velocity- 
gradient tensor dv a /dx/3 in an MD cell but can reproduce 
a velocity profile for an arbitrary (symmetric) strain- 
rate tensor E a p using three components of the velocity- 
gradient tensor, e.g., dv x /dy, dv z /dy, and dv x /dz. For 
simple fluids, the antisymmetric part of the velocity- 
gradient tensor, fl a p = ^(dv a /dxp — dvp/x a ), does not 
affect the stress tensor T a p. Thus, the present multi- 
scale method using the Lees-Edwards boundary condi- 
tion in each MD cell is applicable to the general (three- 
dimensional) flows of simple fluids. 

The off-diagonal stress tensor T' a p is computed accord- 
ing to E' a g and then passed to CFD after transform- 
ing back into the original coordinates, T a p. For two- 
dimensional flows [d/dz=0 and v z =0], and E' are ex- 
pressed as 



8 = 



cos V sin ( 



sin tt cos c 

K y = E' yx = -E xx sin 26 + E xy cos 26, (14) 



(13) 



where 



^itan- 1 



(15) 



Non-equilibrium MD simulations for simple shear flows 
in the rotated Cartesian coordinates are performed in 
many MD cells according to the local strain rate E n s 
defined at each lattice node of the CFD. Once a local 
stress tensor P' a a is obtained at the MD level, the local 
stress at each lattice node P a p in the original coordinate 
system is obtained by combining the pressure p obtained 
a priori by CFD and a tensor T' a a obtained by subtracting 
the isotropic normal stress components from P' a p as 



p = e T [- P i + T']e = - P i + & T T'e, 



(16) 



where I is the unit tensor. 

In the MD simulations, we solve the so-called SLLOD 
equations of motion with the Gaussian iso-kinetic 
thermostat -£t21 



dr j 
dt 



Pi , • 

— +-yr yi e x 

m. J 



dt 



fj - iPv <- 



(17a) 



(17b) 



where e x is the unit vector in the x direction and the 
index j represents the jth particle (j = 1, • • • ,N). Tj 
and Pj + •m'^ry^e. x are the position and momentum of 
the jth particle, respectively, fj is the force acting on 
the jth particle due to the LJ potential in Eq. (j4]), and 
7 is the shear rate experienced by each MD cell, which 
is written as 7 = 2E' xy in the present multiscale scheme. 
Note that, in the SLLOD equations, Pj/m represents the 
deviation of velocity of each particle from the mean flow 
velocity ^r y ^e x in the MD cell. The friction coefficient 
C is determined to satisfy the constant temperature con- 
dition dT/dt = with T = Y,j \pj\ 2 /3mN. The friction 
coefficient £ is calculated as 



(18) 
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We integrate Eq. (JT7|) with Eq. (JTSJl 

using the 

leapfrog algorithm 3 - 3 , with the Lees-Edwards sheared pe- 
riodic boundary condition in the y direction and periodic 
boundary conditions in each cubic MD cell. 

The space integral of the instantaneous microscopic 
stress tensor reads as 

vp af> =-22p-iPfj - E — r (19) 

j=l allpairs 

where we rewrite the momentum of the jth particle, 
Pj + mjr y .e x , as pj. Here, V is the volume of the MD 

cell, i.e., V = Z MD , and £ is the relative vector rj — Ty 
between the two particles r,- and ry . The density is also 
fixed in each MD simulation by the number of particles 
and the box size. Thus, we calculate the local stress using 
so-called NVT-ensembles, P'(p,T,E'). The macroscopic 
stress is averaged in steady states after the transient be- 
haviour has vanished. In the present computations, the 
transient time in each MD process is set as a tenth of the 
CFD time step, At/10. The initial states in each MD 
simulation are created by re-scaling the thermal velocities 
of molecules to be fixed to the local temperatures with- 
out changing the molecular configurations from those ob- 
tained by the previous process. 

In the following, At and Ax represent the time-step 
and the mesh size of CFD calculations, and £md and ^md 
represent the sampling duration and the side length of 
a MD cell, respectively. The two parameters At/tMD 
and Ax/Imd represent the efficiency of our multiscale 
simulations. In the present simulations, we fix the time 
step of the MD simulation A£md as AtMD=0.005. The 
sampling durations of the MD simulation £md are set to 
be larger than the correlation time of the temporal shear 
stress for the bulk fluids, £md ^ tlj- The time step of 
CFD At is set to be small enough for the CFD solutions 
to be stable. 



B. Driven cavity flows 

The simple LJ liquid is contained in a square box with 
a side length L. At t=0, the upper wall starts to move 
from left to right at a velocity V w . A non-slip bound- 
ary condition is applied at each wall: v x —V w and v y —0 
at y=L and v x =v v =0 at the other walls. At the up- 
per left and right corners, v x —V w and v y =0 are applied. 
Although the non-slip condition causes singularities in 
the strain rates and stresses at the upper corners in the 
continuum description, we use the conventional non-slip 
boundary condition even at the corners to compare the 
results to those obtained by the usual CFD calculation. 34 
We note that multiscale simulations to resolve the singu- 
lar forces at the corners can be found in Refs. I35ll36l . 
The multiscale simulations are performed for the param- 
eters listed in Table fl] Here, the Reynolds number is 
defined as pUL/rj, and the viscosity of the corresponding 
LJ liquid is r]=l. 7, which is calculated by the integral 



Imp tup N L U Re Ax /Imp At/tup 



CI 


6.84 


4.68 


256 


219 


0.46 


59 


1.0 


1.0 


C2 


6.84 


9.36 


256 


438 


0.91 


235 


2.0 


2.0 


C3 


6.84 


3.11 


256 


876 


1.83 


941 


4.0 


4.0 


C4 


8.62 


9.36 


512 


438 


0.91 


235 


1.59 


2.0 


C5 


8.62 


18.7 


512 


438 


0.91 


235 


1.59 


1.0 


C6 


10.86 


15.7 


1024 


2780 


0.58 


941 


8.0 


8.0 



TABLE I: Simulation parameters for cavity flows, /md , £md , 
and N are the side length of MD cell, sampling duration of 
MD simulation, and number of particles contained in each 
MD cell, respectively. L, U , and Re are the system size, ve- 
locity of the upper wall, and Reynolds number on the system, 
respectively. Ax /Imp and At/tMD are the efficiency factors 
of the present multiscale simulations. 



of the stress-relaxation function G(t) shown in Fig. [5J 
rj = J Q G(t)dt. The computational domain is divided 
into 32x32 uniform lattices. 

The results of the multiscale simulations are shown in 
Figs. 0] and [5] Figure H] shows the steady-state velocity 
profiles time-averaged over t/At= [900,1000]. It is clear 
that our multiscale method can successfully reproduce 
the characteristic flow properties of cavity flows with 
different Reynolds numbers; the size of the vortex be- 
comes larger as the Reynolds number increases. Figure 
[5] shows the snapshots of velocity profiles for the case 
of Re=980 at different time steps. Here, the results ob- 
tained by multiscale simulations with the efficiency fac- 
tors A:e/Zmd=A£/£md=8 are compared with the results 
obtained for the Newtonian fluid. The results show that 
a small vortex first appears in the upper-right corner 
and then moves gradually toward the centre of the box, 
with the size of the vortex increasing as time passes. Al- 
though fluctuations are seen in the instantaneous velocity 
profiles, our multiscale method can successfully repro- 
duce the characteristic behaviours in the time-evolution 
of driven cavity flow. 

The comparisons of the spatial variations of local ve- 
locities, strain rates, and stresses obtained by the mul- 
tiscale simulation and those of the Newtonian fluid are 
shown in Figs. [Bland [7] It is seen that the deviations of 
the local shear stresses obtained by the multiscale sim- 
ulations are larger than those of the local velocities and 
strain rates. The instantaneous fluctuations of the local 
shear stresses are notable, but they are smoothed by tak- 
ing the time averages. The fractional RMS deviations of 
the local velocities and stresses obtained by the multi- 
scale simulations from those obtained by the usual CFD 
are also shown in Table |H] The comparisons between 
C2, C4, and C5 show the effect of the efficiency factors 
At/tMD and Ax/Imd on the RMS deviations with fixed 
system parameters L, U, and Re. The RMS deviations 
increase as the efficiency factors increase with fixed sys- 
tem parameters, e.g., L, U, and Re. On the other hand, 
in the comparisons between CI - C3, the RMS devia- 
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FIG. 4: The steady-state velocity profiles of the cavity flows for (a) Re=59 (CI in Table |TJ, (b) Re=235 (C2 in Table H]), and 
(c) Re=941 (C3 in Table [[J .The velocity profiles are time averaged over t/Ai=[900,1000]. 



RMS deviation 



CI 
C2 



0.049 
0.034 
C3 0.024 
C4 0.024 
C5 0.017 
C6 0.022 



0.36 
0.25 
0.39 
0.17 
0.12 
0.36 



TABLE II: Fractional root-mean-square (RMS) de- 
viations in the steady states, which are defined as 

Jjf dx 2 J T dt(Q - Q N s) 2 /L 2 r/Max|Q N s|, for the velocity 
(Q=v) and for the shear stress (Q=P xy ). Here, T represents 
the sampling duration in the calculation of the RMS in the 
steady state. 



tions of the velocities slightly decrease in the order CI 
- C3, although the efficiency factors increase in the or- 
der C1-C3. This unexpected finding is explained by the 
fact that the local strain rates increase as the Reynolds 
number becomes larger, and thus the local shear stresses 
are also large in comparison to the noise intensity. The 
comparison of C3 and C6 also shows that the fractional 
RMS deviations decrease even though the efficiency fac- 
tors increase. It should be noted that, in the simulation 
parameter C6, the number of particles contained in each 
MD cell is four times greater than in C3 and the sampling 
duration of the MD simulation is about five times longer 
than in C3. The values of C3 and C6 in Table [TT] can 
be explained by the property of fluctuations discussed 
below. These facts indicate that multiscale simulations 
can be performed with high efficiency factors for large 
Reynolds numbers and large system sizes. 

As mentioned above, the ratios A£/£md and Ax/Imd 
measure the efficiency of our multiscale simulations. For 
example, in the case of A</<md=A:z;/^md=8, the compu- 
tational efficiency is approximately 8 2 x 8 times better 



than that of a full MD simulation. The larger the ratios, 
the more efficient the simulations; however, the statistical 
fluctuations also increase. In the following part, we will 
discuss the nature of the fluctuations in more detai h 37 i 38 
To handle the statistical noise explicitly, we rewrite Eq. 
d]) as 



p = - P i + e T (Ti + R')e, 



(20) 



where the off-diagonal stress tensor T", which is to be de- 
termined by MD sampling, is decomposed into the non- 
fluctuating stress Tl and the fluctuating random stress R! 
due to the thermal noise. The magnitude of each com- 
ponent of the random stress included in MD sampling, 
(^MDpq) where p and q represent the index in Cartesian 
coordinates and do not follow the summation convention, 
should depend both on the size of the MD cell Imd and 
the length of time £md over which the average is taken at 
the MD level; (R' MDpq ) = (R pq (l M D,t M D) 2 ), where R(l,t) 
represents the random stress tensor averaged in a cubic 
with a side length I and over a time duration t. 

At the CFD level, which is discretised with a mesh size 
Aa; and a time-step At, the physically correct magnitude 
should be (R^ FBpq ) = (R pq (Ax, At) 2 ) . If the central-limit 
theorem, (R pq (l,t) 2 )(x l/l 3 t, is assumed, the following 
simple formula can be used: 



Ax 
Imd 



At 



Oml> 



CFDpg/ • 



(21) 



This line of reasoning leads to the following expression 
for the correctly fluctuating stress tensor P: 



P=-pI+ Q 1 




'md \ J / £md \ n r 
Ax") {'AT) KMB 



e (22) 



to be used in CFD instead of Eq. (fl6|) . This equation 
indicates that if we could re- weight the randomly fluctu- 
ating part R' while the non-fluctuating part T£ remains 
untouched, hydrodynamic simulations including correct 
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Hybrid Navicr- Stokes 




(c) t/Ai=600 



FIG. 5: Time evolutions of the velocity profile for the cavity flow with Re=941. The left column shows the results of the 
multiscale simulation with Ai//md = At/tMD = 8 (C6 in Table [I]) and the right column shows the corresponding CFD results. 



thermal fluctuations can be performed within the present 
framework. Incidentally, we can explain the results of the 
fractional RMS deviations in Table [XT] using the relation 
(f2Tj) while noting that the quantity R is normalised by 
P U\ 



transformation of T' xy is defined as 



2M-1 2Af— 1 



K y {k} = ^ ]T T^ y {x}eM-ik-x), (23) 



n x — n y —0 



We note that the important key toward the develop- 
ment of the multiscale simulation is the separation of 
and R' . We thus carried out the spectral analysis for the 
fluctuations in the total stress tensor computed directly 
from MD simulations, T' = T' + R'. The discrete Fourier 



where x = (n x Ax, n y Ax) is the position of each lat- 
tice node (n x ,riy), k — (2itm x /L,2Trm y /L) is the wave 
vector, n x ,n y ,m z ,m y are integers, M is the lattice 
number in each x- and y-axis, and T xy {x} is defined 
as f' xy {x}=T xy {x + Ax/2,y + Ax/2) for < x,y < 
L, f' xy {x}=T' xy {2L - x,y} for L < x < 2L, and 
T'M=T'{x, 2L - y} for L < y < 2L. 
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FIG. 6: Comparisons between the instantaneous profiles of velocity Uv a , strain rate (U/L)E xy , and shear stress {pU 2 )P xy on 
the line of a;/L=0.5 obtained by the multiscale simulation with Ax/Imo = At/tMV = 8 (C6 in Table|T| and those obtained by 
CFD for the cavity flow for i?e=941. The symbols show the results of the multiscale simulation and the solid lines those of 
CFD. The dashed lines in (b) and (c) show the time-averages over t/ At— [900, 1000]. 



1 1 1 1 1 1 1 1 1 1 
Peases V 




- 




Vx 




■ ■■■■■■ 







0.2 0.4 0.6 0.8 1 

x/L 
(a) 

FIG. 7: Comparisons of the instantaneous profiles of velocity Uv a , strain rate (U/L)E xy , and shear stress (pU 2 )P xy on the 
line of y/L=0.5 obtained by the multiscale simulation with Ax/Imd = Ai/tMD = 8 (C6 in Table|I]) and those obtained by CFD 
for the cavity flow for J?e=941. See also the caption in Fig. [6] 



The power spectra {\n' xy {k}\ 2 ) calculated from our 
multiscale simulations are plotted in Fig. [S] (a) for the 
case of Ai/iMD = Ai/!md = 1, which corresponds to 
the case of Fig. [2(a) (CI in Table 0}. The angle bracket 
(• • • ) indicates the time average taken at steady state at 
the CFD level. One can see that the overall structure is 
rather simple. There exist a relatively large peak around 
k = and rather flat distributions throughout the k 
plane. The former corresponds to the contributions from 
the non- fluctuating part Tl, and the latter corresponds 
to the contributions from the random stress R' . The 
same quantity obtained by conventional fluctuating hy- 
drodynamics using a constant Newtonian viscosity and 
a random stress, the intensity of which is determined by 
the fluctuation-dissipation theorem^ 7 , is shown in Fig. [8] 
(b) for comparison. 3 - Those two plots are surprisingly 
similar , including the fluctuation part. This similarity 
means that our multiscale simulation generates fluctu- 
ations quite consistent with fluctuating hydrodynamics 



using the fluctuation-dissipation theorem in the case of 

Az/Z M D=At/*MD=l. 

In Fig. UU one sees how the fluctuations depend on 
the ratios Ax/Imd and At/twin- Here, in comparison to 
the reference case (a) [Ax/ZMD=At/iMD=2], the num- 
ber of particles used in the MD simulations are dou- 
bled in the case of (b) [Ax/Imd=1-59, A£/iMD=2], and 
both the number of particles and the sampling duration 
used to take the time average are doubled in the case 
of (c) [Az/Zmd=1-59, Ai/iMD=l]- The noise intensity 
decreases with decreasing ratios Ax/Zmd and At/tMU, 
consistent with the central-limit theorem Eq. (f2"Tj) ; i.e., 
the noise intensity in (b) is approximately half that in 
(a), and the intensity in (c) is about one quarter of that 
in (a). 

From the overall properties, we can confirm that the 
fluctuations arising in the present multiscale simulations 
are white noise, which obeys the central-limit theorem, 
written as Eq. ([2~Tj) or (|2"12j) . The results also suggest that 
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FIG. 8: The fluctuations of T' xy for the case of cavity flow with Re=59. The power spectrum (\II' xy {k}\ 2 } for the present 
multi-scale model with Ax/Imu = At /tun = 1 is shown in (a); the corresponding result from the fluctuating hydrodynamics 
is shown in (b) for comparison. U' xy represents the discrete Fourier transform of T' xy . m a is defined as m Q =(L/27r)fc a , where 
k is the wave vector. The insets on each figure show the (\n' xy \ 2 )-m x plane. 




(a) 



(b) 



(c) 



FIG. 9: The fluctuations of T xy for the case of cavity flow with Re=235. The power spectra (|77^ a {fe}| 2 ) is plotted in (a) for 
the case of Fig. [4] (b) . In (b), only the number of particles is doubled; other parameters are unchanged from (a). In (c), both 
the number of particles and the sampling time of T' xy are doubled. U ~' xy represents the discrete Fourier transform of T xy . m a 
is defined as m a — (L/2Tr)k a , where A; is the wave vector. The insets on each figure show the (\n' xy \ 2 )-m x plane. 



some white-noise-reduction algorithm, such as a low-pass 
filter, could be useful in the CFD-MD coupling processes. 



III. ONE-DIMENSIONAL POLYMERIC FLOWS 

As we have seen, the relaxation time of the local 
stresses in simple fluids is very short; for the LJ liquid in 
the previous section, it is estimated as tlj < 0.1 in the 
LJ unit time (Fig. [2]). Thus, we can predict the local 
equilibrium state at each instant in the sampling of the 
local stresses for the simple fluid (Fig. [3]) . For polymeric 
liquids, however, the relaxation time is much longer than 
that of simple fluids, and thus it happens to be larger 
than the characteristic time of the macroscopic flows. In 
Fig. [3J we show the stress relaxation of a model poly- 
mer melt, which is discussed in this section. In this case, 
we cannot predict the local equilibrium states within the 
time-step duration to resolve the macroscopic flow be- 
haviours. Instead, we must consider the history of the 
slow dynamics of polymer chains in each fluid element, 



i.e., the "memory effect". 

In this section, we consider the flow behaviours of a 
polymeric liquid confined in parallel plates; we assume 
that the flow direction is restricted to being parallel to 
the upper and lower plates and that the local macro- 
scopic quantities, e.g., the flow velocities, strain rates, 
and stresses, are uniform along the flow direction. This 
assumption allows us to calculate the macroscopic quan- 
tities using the usual fixed-mesh system without tracing 
the advection of a fluid element that contains a mem- 
ory of the configuration of polymer chains. We note 
that one must treat the convection of the memory along 
the stream lines in the general two- or three-dimensional 
flows. The extension to polymeric flows with the advec- 
tion of the memory will be given in the next section. 

A. Multiscale modelling and the simulation 
method 

We consider a polymer melt with uniform density po 
and temperature Tq between two parallel plates (Fig. 
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FIG. 10: Schematics for the geometry of the problem, mesh system, and time-evolution scheme for one-dimensional flows of 
polymeric liquids in parallel plates. 



[TUJa)). Both plates can move in the x-direction. The 
melt is composed of short chains with ten beads. The 
number of beads comprising each chain in the MD sim- 
ulation is represented by Nb- Thus, Nb = 10. All of the 
beads interact with a truncated Lennard- Jones potential 
defined by Eq. ((4]) in Sec. By using only the re- 

pulsive part of the Lennard- Jones potential, we prevent 
the spatial overlap of the particles. Consecutive beads 
on each chain are connected by an anharmonic spring 
potential 



U F (r) 



; k c R 2 \n [l-(r/R ) 2 



(24) 



where k c =30e/a 2 and Ro=1.5a. The number density of 
the beads is fixed at po / m=l j <r 3 , where m is the mass of 
the bead particle. With this number density, the con- 
figuration of bead particles becomes severely jammed 
at low temperatures, resulting in the complicated non- 
Newtonian viscosity and long-time relaxation phenomena 
typically seen in glassy polymers 

We assume that the macroscopic quantities are uni- 
form in the x- and z-directions, d/dx—d/dz=0. The 
macroscopic velocity v a is described by the following 
equations: 



Po- 



dv x dP, 



X1J 



dt 



dy 



(25) 



and v y —v z —0, where t is the time and P a p is the stress 
tensor. Here and afterwards, the subscripts a, j3, and 
7 represent the index in Cartesian coordinates, i.e., 
{a,(3,j}={x,y,z}. We also assume a non-slip boundary 
condition on each plate. 

We use a usual finite volume method with a staggered 
arrangement for the CFD calculation; in this method, 
the velocity is computed at the node of each slit and 
the stress is computed at the centre of each slit. For the 
time-integration scheme, we use the simple explicit Eulcr 
method with a small time step At. A local stress is calcu- 
lated at each time step of CFD using the non-equilibrium 
MD simulation with a small cubic MD cell with a side 
length Imd associated with each slit according to a local 



strain rate. In each MD simulation, we solve the SLLOD 
equations of motion with the Gaussian iso-kinetic ther- 
mostat, Eq. (fTTf with (fT5|) . as in the previous section. 
The space integral of the microscopic stress tensor reads 
as 



k=i 3 =i 

EE 

fc=ii=i 



allpairs 



(26) 



where we rewrite the momentum of the jth bead on the 
kth chain, p k +m A fr y k e x , as p*. Here, the indexes k and 
j represent the fcth polymer chain (k = 1, • ■ • , N p ) and 
the jth bead (J = 1, • • ■ , Nb) on each chain, respectively. 
N p and iVf, represent the number of polymer chains and 
beads on each chain, respectively. £ in the right-hand side 
ofEq. (|2"6"]l represents the relative vector rj* — r j, between 

the two beads, r k and r|, , in the second term and the 



relative vector r. 



'i+i 



between the two consecutive 



beads on the same chain, r| and in the third term. 

In the present problem, we cannot assume a local equi- 
librium state at each time step of the CFD simulation be- 
cause the relaxation time of the stress may become much 
longer than the time step of the CFD simulation (with 
which the macroscopic motions should be resolved). In 
the current simulations, the simple time averages of tem- 
poral stresses of the MD (averaged over the duration of a 
time-step of the CFD simulation) are used as the stresses 
at each time step of the CFD calculation without ignor- 
ing the transient time necessary for the MD system to be 
in the steady state (Fig. [TU](c)). Thus, the time integra- 
tion of the macroscopic local stresses P a p is performed 
with the instantaneous microscopic stress tensor cr X y(t) 
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as 



t+At 



/ 3 
MD Jt 



P a p(t',y)dt' 

t+At 

CTQ^(r;7(t,y))dr, 



(27) 



where At is the time step of the CFD simulation and 
!md is the side length of a cubic MD cell. The second 
argument of o a p in Eq. ([27]) . i.e., j{t,y), is constant in 
the integral interval, which indicates that the shear rate 
to which each MD cell is subjected is also constant over 
a duration At in each MD simulation. The final configu- 
ration of the molecules obtained at each MD cell at each 
time step of the CFD calculation is retained as the ini- 
tial configuration for the MD cell at the next time step of 
the CFD. Thus, we trace all of the temporal evolutions 
of the microscopic configurations at the MD level so that 
the memory effects can be reproduced correctly. We can 
reduce the computation time needed for the spatial in- 
tegration compared to that in a full MD simulation by 
using MD cells that are smaller than the slit size used in 
the CFD simulation. The performance efficiency of the 
present multiscale simulation is represented by a saving 
factor defined as the ratio of the slit size used in the CFD 
simulation Ax to the cell size used in the MD simulation 
Imd, Ax/Imd- It also should be noted that, in addition to 
the saving factor Ax/Imd, the present multiscale method 
is quite suitable for a parallel computational algorithm 
because the MD simulations associated with each mesh 
of the CFD, which consume a large portion of the total 
simulation time, are performed independently. 

We measure the space and time in the units of a and 
To = yj m<r 2 /e, respectively, just as in the previous sec- 
tion. The temperature T is measured in the unit of e/ks- 
In the following simulations, the temperature T and den- 
sity p of the melt are uniform and fixed as T = 0.2 
and p — 1.0. At this number density and low tem- 
perature, the polymer melt involves complicated non- 
Newtonian rheology. The time step of the CFD simu- 
lation At, the sampling duration of the MD simulation 
iMD, and the time step of the MD simulation At are 
fixed as Ai=t M D=l and Ar=0.001. Thus, 1000 MD 
steps (M = 1000 in Fig. [KjJ) are performed in each 
MD cell at each time step of the CFD computation. One 



hundred chains with ten beads are confined in each cubic 
MD cell with a side length ^md=10; thus, iV p =100 and 
N h =l0. 



B. Oscillatory flows 

The lower and upper plates start to oscillate with a 
constant angular frequency ujq at a time t = as, respec- 
tively, 



v w = ±w cos(w t), 



(28) 



with vq = TqljqH. Here, T represents the amplitude of 
strain on the system. In the following simulations, the 
width between the plates, 2H, is fixed as 277=5000. In 
the CFD calculations, the lower half of the volume be- 
tween the plates is divided into 128 (or 64) mesh intervals 
with a mesh size of Ax=19.5 (or 39.0), and the symmetric 
condition is used at the middle between the plates. The 
128-mesh system is used for a large To, and the 64- mesh 
system is used for a small Tq. Thus, the saving factor 
Ax/Imd in the multiscale simulations is Ax/7md=1-95 
(or 3.9). 

When the velocity amplitude of the oscillating plates 
vo is small enough, ti «l (or T -C 1), one can assume 
the linear response of local stress to the local strain rate 
as 



P, 



G(t - t')j(t')dt', 



(29) 



where G{t) is the stress-relaxation function, which is 
shown in Fig. [2j By taking the Fourier transform of 
Eq. with Eq. ([2ljp, we obtain 



dy 2 



ip uj 
1 

r/* 



(30) 



where is the Fourier transform of the local velocity, 
v x(y) = I^°oo v ft, y)e~ iult dt, and rj* is the complex viscos- 
ity, rj* = if + vq" , obtained by the Fourier transform of 
the stress-relaxation function as r]*(uj) = f Q °° G(t)e~ lut dt 



Together with the boundary condition, Eq. (|28|) , the an- 
alytical solution for the linear response regime is 



J 



v x (y,t) = v {e~ k °y eMK^ot ~ k'Jy)} _ g-^-i/) ex p[i( Wo t - k$(2H - y))}} (l 



-2kXH 



(31) 



with 



K 



(32a) 



and ko = k' Q + iky, where 9 is the loss tangent defined 
as 6 = tan _1 (ry"/ry'). The complex viscosities rj and rj" 
can be calculated from the numerical data of the stress- 
relaxation function G(t) shown in Fig. [5] In the present 
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paper, the analytical solution of Eq. (f3Tj) with Eq. (|32|) is 
used to show the validation of our multiscale simulations. 
We will demonstrate that the results of the multiscale 
simulation become closer to the analytical solutions as 
the strain amplitude Tq decreases. 

Figures [TT] and Q2] show the velocity profiles of the 
polymer melt in steady oscillations. It is seen that, as 
the amplitude of the oscillating plates Tq increases, the 
thickness of the boundary layer becomes thinner. This 
phenomenon is related to the shear thinning occurring 
near the oscillating plate, where the local strains are large 
enough for the local viscosities to deviate from those in 
the linear-response regime. 

With use of the Fourier transform for the time evo- 
lutions of local velocities and selecting the mode of the 
oscillation frequency of the plate ojq, we can express the 
time evolution of the oscillation velocity as 

v x {y,t) = \v x \(y)cos(ut + <f>(y)), (33) 

Here, the contributions of odd higher harmonics, 3ujq, 
5u!q, and • • • , are ignored because they are much smaller 
than the mode of wo- 15 A detailed investigation of the 
nonlinear rheology of the present melt can be found in 
Ref. 0. Figure [T3] shows the comparison of the ampli- 
tudes of local velocity \v x \ between the various ampli- 
tudes of strain on the system Tq. As the amplitude of 
oscillating plates Tq decreases, the profiles of the ampli- 
tudes of the local velocity 1^1 become closer to that for 
the linear response. At large values of To, |ttc| rapidly 
decreases near the oscillating plate and approaches the 
linear-response profile as the distance from the plate in- 
creases. One can also measure the wavelength for the 
velocity profile, A, using the relation 0(A) = 0, where 
(j>(y) is the local phase retardation in Eq. (|3"3")l . For the 
analytical solution Eq. (|3T|) . the wavelength A is written 
as A = 2-K/k.Q. Figure [TH shows the comparisons of the 
wavelengths A between various strain amplitudes for the 
system Tq. It is seen that the wavelengths are very close 
to those for the linear response at all oscillation frequen- 
cies wo at r =0.005 but that they deviate more strongly 
from those for the linear-response regime as Tq increases. 
The deviation is slightly larger because the oscillation 
frequency loq is larger at the same value of Tq. 

We also measure the "local" viscoelastic properties in 
terms of the "local" complex viscosity if (=77' + if]"). The 
real and imaginary parts of the complex viscosity, rf and 
77", represent the viscous and elastic responses, respec- 
tively, in the shear stress. The complex viscosities rj' 
and 77" can be calculated as follows. As in the case of 
the oscillation velocity, by using the Fourier transform 
of the time evolution of the local strain rates and se- 
lecting the mode of the oscillation frequency loq, we can 
reconstruct the time evolution of the local strain rate 
7 in the form of 7 — \j\(y) cos(woi + ip(y))- In the 
same way, we can also write the local shear stresses as 
p xy = Kyiv) cos(ui t<ip) - P" y (y)cos(w o t-0)- The com- 
plex viscosities rj' and 77" are obtained by rj = Pxy/\j\ 
and 77" = P^ y /I'll' respectively^ 



Figure [15] shows the spatial variations of the local com- 
plex viscosities 77' and 77" and the amplitude of the lo- 
cal strain |7| (=|7|/wo) for the oscillation frequencies 
ljq = 6.1 x 1CP 3 and 0.025 for various strain amplitudes 
on the system Tq. For comparison, the results for the 
linear response are also shown. As the distance from 
the plate increases, the amplitude of the local strain | — / 1 
decreases monotonically. When I7I is smaller than a 
few percent, the local complex viscosities become con- 
stant, and these constant values correspond to those for 
the linear response. As the amplitude of the oscillat- 
ing plates decreases, the profile of the amplitude of local 
strain I7I approaches that for the linear response. In 
addition, the range in which the local complex viscosi- 
ties agree with those for the linear response broadens 
near the plate. In particular, the dynamic viscosity 77' in 
Fig. I15f a) is almost uniform over the entire domain at 
Tq = 0.005. Shear-thinning is observed in the vicinity 
of the oscillating plate. That is, the complex viscosities 
77' and 77" decrease near the oscillating plate, where the 
local strain rapidly increases. The elasticity 77" decreases 
more rapidly near the oscillating plate than the viscosity 
77'; thus, the viscosity 77' is dominant in close vicinity to 
the plate, while the elasticity 77" is larger than the viscos- 
ity 77' far from the plate. In the case of a high oscillation 
frequency and large amplitude of oscillating plates, say 
cjo=0.025 and Tq = 0.5, the melt behaves as a viscous 
liquid in the vicinity of the plate because the viscosity 
rj' is much larger than 77", 77' 3> 77", but as a viscoelastic 
solid far from the plate because the elasticity 77" is larger 
than 77', 7/ > rj". In the transient region, the melt be- 
haves as a viscoelastic liquid. Thus, the melt has three 
different rheological regimes of viscous liquid, viscoelastic 
liquid, and viscoelastic solid over the oscillating plate at 
high oscillation frequencies and large amplitudes of the 
oscillating plate. 

The crossover of the local rheological properties can be 
characterised by the local Deborah numbers De defined 
by the oscillation frequency of the plate loq and the local 
relaxation times of the monomer structures and polymer- 
chain configurations, r Q and t r (which are referred 
to as the a relaxation time and the Rouse relaxation 
time, respectively), as De a = luqt and De R = ^>qt R - 
The relaxation times r Q and t r decrease as the strain 
rate 7 increases. Thus, the local Deborah numbers are 
smaller closer to the oscillating plate. We found that the 
crossover of 77' and 77" occurs at the position where the 
local Deborah number De a is unity. In the viscous liquid 
regime, the local Deborah number De R is less than unity. 
For more details, consult Refs. [TillTsl . 



IV. ADVECTION OF STRAIN-RATE MEMORY 
IN TWO-DIMENSIONAL POLYMERIC FLOW 

In the previous section, we described "the memory ef- 
fect" in one-dimensional polymeric flows. In general two- 
dimensional polymeric flows, we must consider the advec- 



14 




v x /v v x /v v x /v 

(a)r ^o (b) r = 0.005 (c)r = o.i 

FIG. 11: Velocity profiles of the melt at each 7r/4 phase shift for wo = 6.1 x 10~ 3 . (a) The analytical solution with the 
numerical data of rj' and rj" calculated by the Fourier transform of G(t) shown in Fig. [2] (b) and (c) The results of the 
multiscale simulations for ro=0.005 and 0.1, respectively, where instantaneous velocity profiles at each fixed phase are averaged 
in the steady oscillation. 
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FIG. 12: Velocity profiles of the melt at each 7r/4 phase shift for uio = 0.025. See also the caption in Fig. QT] 
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FIG. 13: Comparison of the amplitudes of local velocity 1^ | for various strain amplitudes on the system To. (a) ujo = 6.1 x 10 3 
and (b) uj = 0.025. 
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tion of the strain-rate memory that is preserved in mi- 
croscopic degrees of freedom in a fluid element. Because 
the general two-dimensional flow is inhomogeneous along 
the flow direction, the advection of the memory affects its 
flow behaviour when the characteristic relaxation time r 
of a polymeric liquid is larger than the transit time Tt 
that the polymeric liquid takes to move through a dis- 
tance equal to the size of a fluid clement. Polymers in 
the fluid element experience a local strain rate during 
T t . When a strain rate is imposed on the fluid element 
during r t (< r) in the inhomogeneous flow, the stress 
response of the fluid element is different from that in 
the steady state under the same strain rate because the 
strain-rate memory in the stress persists during r longer 
than Tt. Therefore, the stress does not reach the value 
at the steady state. Because the transit time becomes 
infinitesimal as we decrease the fluid-element size to in- 
crease the spatial resolution of the numerical simulation, 
consideration of the advection of the memory is essen- 
tial for a polymeric liquid with a finite relaxation time in 
inhomogeneous flow fields. 

To deal with the advection of the memory in a manner 
consistent with the macroscopic flow, we employ a La- 
grangian fluid method in which microscopic simulators 
are assigned to the fluid elements^ Because each fluid 
element moves along the flow direction, the microscopic 
variables are advected in a Lagrangian manner. 

In this section, we discuss the flow of a well-entangled 
polymer melt that is inhomogeneous in the flow direc- 
tion. The polymer melt discussed in the previous sec- 
tion was composed of short polymer chains without en- 
tanglements among polymers. In the short-chain sys- 
tem, the memory effect is mainly caused by the relax- 
ation of the polymer orientation. In the polymer melt 
flow without entanglements, the relaxation time of the 
polymer orientation can be considered to be the same as 
the relaxation time of the polymer's extension because 
the origins of these relaxations are same. That is, both 
relaxation times are determined by the competition be- 



tween the following two forces: i) an elastic force coming 
from the chain-conformational entropy and ii) a frictional 
force coming from the interaction between beads. In an 
entangled polymer melt, however, these two character- 
istic relaxations are fundamentally different because of 
the presence of the entanglements. Therefore, in deal- 
ing with entangled polymer melts, we will observe the 
effects of each of these different relaxations on the flow 
behaviour. We will ultimately conclude that the memory 
effect relates not only to the relaxation of orientation but 
also to the relaxation of extension^ 



A. Lagrangian fluid model with microscopic 
internal degrees of freedom 

We consider a polymer melt composed of well- 
entangled linear polymers. In this polymer melt, the 
number of beads between two adjacent entanglements 
along the polymer chain has been found to be approxi- 
mately 35 in terms of the MD polymer chain^S. Therefore, 
to describe the entangled polymer melt, each polymer 
chain in the MD simulation must consist of more than 100 
beads. Dealing with the entangled polymer melt in the 
multiscale simulation with MD simulators seems impos- 
sible because of the extremely high numerical expense. 
To decrease the numerical cost of the multiscale simu- 
lation, we employ a coarse-grained slip-link model that 
reduces the 35 beads to a single primitive path. Specif- 
ically, by using the coarse-grained model, the degrees of 
freedom become less than l/35 th -of those in the corre- 
sponding MD simulation. The great reduction in the nu- 
merical expense using the coarse-grained model enables 
us to simulate the entangled polymer melt in the multi- 
scale approach with almost the same numerical expense 
as in the case of the unentangled polymer melt with the 
MD simulation. 

As in the previous sections, we solve the macroscopic 
momentum equation of fluid 

Pi ^r = v ' °* ~ Vpt + pb ' (34) 

where the subscript i represents the i-th fluid element 
that is handled as a particle at the position r^, pt is the 
density, Vj is the velocity vector, er.j is the stress tensor, 
Pi is the isotropic pressure, and F b is a body force. More- 
over, we need to solve the time evolution of the position 
r, that follows 



to consider the fluid in a Lagrangian manner. Once we 
interpret the momentum equation of the fluid (Eq. (|34p ) 

as 
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FIG. 15: Spatial variations of the local complex viscosities rj (solid line) and rj' (dashed line) and the amplitude of local strain 
| — y | (long-dashed line) for the oscillation frequencies ujo = 6.1 x 10~ 3 (a) and 0.025 (b). The blue, green, and red colours show 
the results for Fo=0.005, 0.02, and 0.5, respectively, in both figures. The black lines are the results for the linear response for 
infinitesimally small strains. 
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FIG. 16: Interaction circles with a cut-off length 2h in two- 
dimensional SPH. The small filled and vacant circles (•, o) 
represent fluid particles. The large solid circle represents the 
interaction range of the i-th fluid particle. 

where Fj is the right-hand side of Eq. (|34l) , time-integral 
schemes developed in MD simulations are available. Here 
we implement the velocity- Verlet algorithm to update 
{rj and {vj. 

To solve Eq. we need to know the derivatives of 

the stress and pressure at the position r,-, (V • <x)j and 
(Vp)j. Because we now know that the stress and pres- 
sure have values only at the positions {r^}, we cannot 
obtain their derivatives in a straightforward manner. In- 
stead of directly obtaining the derivatives (V • <x)j and 
(Vp)i, solving a linear matrix equation composed of in- 
terpolated values of <7i and pi enables us to obtain their 
derivatives. The procedures are explained in the follow- 
ing paragraphs. 

In the Lagrangian fluid simulation, as the positions of 
the fluid particles {r.;} change over time, interpolation 
techniques of physical variables at a certain position be- 
tween the fluid particles are required to obtain the spatial 
derivatives V/ or V-f where / (f) indicates an arbitrary 
scalar variable (vector or tensor). In the Eulerian fluid 



simulation, as the interpolation scheme has been well es- 
tablished, we do not need to consider the interpolation 
scheme. In our Lagrangian fluid dynamics simulation, ac- 
cording to the smoothed particle hydrodynamics (SPH) 
simulation^ we employ the Gaussian kernel interpola- 
tion; 

f(T i )=a d ^ f f(r j )W(\T i -r j \,h), (36) 

where /(rj) represents an arbitrary physical variable at 
the position rj in the fluid, /(rj) is the interpolated vari- 
able of /, a is the unit length in the fluid particle simula- 
tion, and d is the dimensionality of the space. The win- 
dow function W(x, h) is assumed to be a Gaussian shaped 
function of distance x and half- value width h = 1.2a, and 
it follows 

W{x,h) = lvfc( e ~*- e -*) (N ^ } ' (37) 
[0 (otherwise), 

where the cut-off radius is 2h. The normalisation param- 
eter A d equals 1.04823, 1.10081 and 1.18516 for d = 1, 2 
and 3, respectively— To obtain an interpolated variable 
/ of a certain variable / at the position of the i-th fluid 
particle, fluid particles in the circular (two-dimensional)/ 
spherical (three-dimensional) region fij with radius 2h in 
which the i-th fluid particle is placed at the centre are 
accounted for, as shown in Fig. 1161 Hereafter we will 
use variables with tildes to represent interpolated values 
using Eq. ([55]). 

To obtain the spatial derivatives of /(rj), i.e., d a f (a — 
{x, y, z}) at the position rj, we solve the following linear 
equation of the matrix form (see Refs. IlBiZBE^): 
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where C 3,7 = 1 — \5p tl and Jg )7 is Kronecker's delta. 
Here we used a reduced representation; Greek indexes 
represent the coordinate axes and the number of rows 
in Eq. (|38l) is seven in two dimensions and 13 in three 
dimensions. The derivatives of the interpolated variables 
in the left-hand side in Eq. ([38]) are defined as 

dj{vi) = a d Y, f{*j)d*W{\Ti - v^h), (39) 

d a dpf(Ti) ^a d J2 ffa)d a d P WQTi ~ Tj\,h), (40) 

and the interpolated values of positional difference be- 
tween the i-th and j-th fluid particles R = Vi — Yj and 
dyadic RR are, respectively, 

R-=a d Y, Ra W{\T i -T j \,h), (41) 

R^Rfi =a d R a R P W(\r i -r j \,h). (42) 

All values in Eq. (13"8"1) . except for d a f and dpd^f, are 
calculated using Eq. <j36j), and Eqs. (J39j)d421). Solv- 
ing Eq. (1551) using the usual linear algebra, we obtain 
d a f (ct = {x, y, z}). The square matrix in the right-hand 
side of Eq. ()38j) is common among the macroscopic vari- 
ables in the same configuration {r.;}. Therefore, once the 
square matrix is decomposed using LU decomposition, 
we can reuse the decomposed matrix, decreasing the nu- 
merical cost of solving the linear equations (Eq. (|3"81 ) for 
variables in the same fluid particle configuration. 

In Eqs. (|3"4"f and (1331) . the value of a physical field at 
the position of a fluid particle is given by the physical 
quantity that the fluid particle holds, and the deriva- 
tives of the physical field can be obtained through Eq. 
(f38|). The density field, however, is not a physical vari- 
able that a fluid particle possesses but a physical variable 
determined from a configuration of fluid particles. Thus, 
we use the following definition for the density field in the 
same manner as SPH^ 5 - 

Pi = Y J ^iW{\T i -T j \,h), (43) 

where m, is the mass of the i-th fluid particle. Here, 
we impose incompressibility to the polymer melt. To 
maintain the initial uniform density at each fluid particle, 
we assume the pressure pi to be proportional to the local 
density pi£& 

Pi = c 2 p t , (44) 



where c is the sound velocity. When the fluid particles 
fill all the space, repulsive forces act between particles 
through Eqs. ^ and flU}. When we set c = 1.0[a/i ], 
where to is the time unit in the Lagrangian fluid dynamics 
simulation, the fluid is almost incompressible when Re < 
1.0, and thus we use this value. 

As mentioned above, the stress tensor <x of the entan- 
gled polymer melt is described by the dynamic response 
of the microscopic state of the polymers. To efficiently 
obtain the stress tensor a of the entangled polymer melt, 
we employ a coarse-grained model for the polymers, i.e., 
the slip-link model 2,3,5 that can accurately describe the 
dynamics of entangled polymers. In a well-entangled 
polymer melt, a polymer chain can be considered to be 
in a tube, composed of surrounding polymer chains^ In 
the slip-link model, the slip-links substitute for the tube 
that constrains the polymer from moving perpendicular 
to the polymer contour direction but allows the poly- 
mer to move parallel to the contour direction, i.e., the 
reptation motion42 Instead of the tube description, rep- 
resentative constraints, i.e., slip-links, are traced in the 
slip-link model. 

Focusing on a section of the polymer chain between 
two adjacent slip-links along a polymer, both ends of the 
section are fixed by the slip-links, and the chain confor- 
mation of the section has a complex curve. Connecting 
the slip-links creates a primitive path with a complex 
curve. Therefore, the orientation of the part of the poly- 
mer chain between two adjacent slip- links is considered 
to be composed of the following: i) the orientation of the 
primitive path that directly connects the adjacent slip- 
links (i.e., the average orientation of the section of the 
polymer chain in the scale of the primitive path length) 
and ii) a random segment orientation on a much shorter 
length scale than that of the primitive path. This random 
segment orientation can be regarded as isotropic because 
the dynamics of segments smaller than the primitive path 
are much faster than the dynamics of the primitive path. 
The stress tensor a depends on the conformations of 
the microscopic polymers, and the stress tensor cr can 
therefore be divided into the following: i) <r p from the 
primitive paths and ii) <r d from the isotropic part. The 
latter stress tensor is assumed to be er d = ?7dD, where 
D = Vv + (Vv) T because of its isotropy. 

At each time step of the slip-link simulation, the con- 
formations of the slip-links are updated by the following 
four steps: 

1 . Affine deformation according to the macroscopic lo- 
cal velocity gradient tensor n = (Vv) T ; 



18 



2. Change in the contour length of the primitive path 
due to the random motion of the slip-links; 

3. Reptation motion along the contour due to the ran- 
dom motion of the centre of mass of the polymer; 

4. Constraint renewal because a slip-link is removed 
when the slip-link has slipped off from the poly- 
mer chain or created when the one of the ends of 
the polymer chain entangles with the other polymer 
chain. 

For more details, see Refs. For a given chain config- 
uration, the stress tensor is calculated by 



'a/3 



E 



r ia r il3 



(45) 



where a s is the unit length in the slip-link model and r a ia 
is the a-component of the z-th segment vector connect- 
ing adjacent slip-links along a polymer. The number of 
slip-links on a chain is represented by Z . The unit of 
stress cr e in the slip-link model is related to the plateau 
modulus Gn by a c = (15/4)Gn£ The slip-link model has 
two characteristic timescales: the Rouse relaxation time 
tr and the longest relaxation time Td- These character- 
istic times are related to Z as follows: tr = Z 2 t c and 
Td oc Z 3A t c , where t e is the time unit in the slip-link 
simulation^ The contour-length relaxation of a confin- 
ing tube, i.e., the extensional relaxation, occurs on the 
timescale of tr, while the orientation relaxation occurs 
on the timescale of Td . 

Each polymer simulator describing a fluid particle com- 
putes the polymer configurations at each time step, and 
the recorded configurations are used as the initial states 
in the next time step. Typically, the macroscopic time 
unit tmacro and microscopic time unit i m i cr o have a large 
timescale gap, and the macroscopic time unit t macro must 
therefore be divided into Nt m i cro . Because the slip- 
link model used here is sufficiently coarse-grained and 
the time unit t e can be the same as the timescale of 
the macroscopic fluid t maclo = to, we employ f macro = 

^micro = tp.- 



B. Two-dimensional polymeric flow with memory 
advection 



To demonstrate the efficiency of the presented multi- 
scale simulation, we consider a two-dimensional polymer 
melt system in which the flow history can affect the flow 
behaviour. Figure [17] shows the system in which a circu- 
lar object with radius r c = 3a is fixed on the centre of 
the square system with sides 2H, where H is set to be 
H = 15a. We imposed the stick boundary condition on 
the velocity at the perimeter of the circular object and 
periodic boundary conditions on all the physical variables 
used at the boundaries of the system. Around the circu- 
lar object, the polymers in a fluid particle experience a 




FIG. 17: Two dimensional polymeric flow system with a cir- 
cular obstacle. The circular object with a radius r c is fixed 
on the centre of a square system with sides 2H. To describe 
the circular object, fixed particles are evenly placed on the 
circle. Mobile fluid particles flow outside of the circle. The 
inside of the circle is vacant. Periodic boundary conditions 
are imposed at the sides of the square system. 



strain rate that depends on the position along the stream 
line, and the conformations of the polymers turn out to 
be different between the upstream and downstream re- 
gions. The difference in the conformations between these 
two regions can be observed macroscopically in the dis- 
tribution of the stress through Eq. (1431) . 

The polymer melt is composed of monodisperse linear 
polymers with an average number of entanglements on 
each polymer of Z = 7. In the bulk of this system, 
the zero-shear viscosity 77Q and the longest relaxation 
time Td are estimated from the slip-link simulation to 
be t]q ~ 17.5cr c i c and Td = 20(M C - 17 Approximately 900 
fluid particles are evenly placed in the system. The cir- 
cular object fixed in the space is expressed by 24 fluid 
particles placed at the perimeter of the circle; the inside 
of the circle is vacant. Each fluid particle consists of 1000 
polymers, enough to describe the bulk rheological prop- 
erties of a single fluid particle^ Under the body force 
F b /(?7 /a< c ) = (5.0 x 10" 4 , 0), the flow becomes a steady- 
state in about 500t c . The average steady-state flow ve- 
locity in this system is nearly equal to (0.055, 0)[a/t e ] for 
a polymer melt flow with Z = 7 and Tjd/rj = 0.1 and the 
Reynolds number Re = pUr c /rf is 0.2. Under these con- 
ditions, the flow is almost laminar, i.e., the flow between 
the upstream and downstream regions can be symmet- 
ric at the steady-state of low Reynolds number flow in 
the Newtonian fluid. In the regions where D xy > l/rj, 
the Weissenberg number We is larger than 1, and the 
flow shows shear-thinning. The shear-thinning behaviour 
is caused by the anisotropic distribution of the polymer 
chain orientation, i.e., the memory of the strain rate that 
the polymer chains have experienced. The orientation 
relaxation is not accomplished in a flow with a strain 
rate larger than 1/rd within a time interval less than Td . 
Similar to the orientation relaxation, the extensional re- 
laxation is not accomplished in a flow with a strain-rate 
flow larger than 1/tr within a time less than tr. Be- 
cause the Rouse relaxation time tr is much smaller than 
the longest relaxation time Td, the extensional relaxation 
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progresses easily compared to the orientation relaxation, 
which is obstructed by entanglements. Therefore, the 
extensional relaxation is difficult to observe in the poly- 
mer melt flow without careful treatment and observa- 
tion. In the present multiscale simulation following the 
Lagrangian fluid dynamics at the macroscopic level, we 
can consider not only the orientation relaxation but also 
the extensional relaxation caused by the memory of the 
strain rate because the microscopic simulators are com- 
posed of entangled polymers and the advection of the 
entire polymer configurations is involved. 

Figures [T5] and [T5] show the transient flow behaviour 
of magnitudes of (a) velocity |v|[a/f c ], (b) strain rate 
\D xy \[l/t c ], and (c) shear stress over zero-shear viscos- 
ity | <r X y | /r]°[l /t c ] at t = 100,200,400, and 800[t o ] in the 
polymer melt and the Newtonian fluid, respectively. To 
obtain these contour maps, we performed a linear inter- 
polation to transform the data at the particle positions 
into values at regular lattice points. To decrease the spiky 
noise of the data, median-filtering, replacing the value 
with the median value obtained from the values of the 
neighbouring pixels and the pixel's own value, was carried 
out to the evaluated data on the regular lattice points. 
Comparing Figs. [18] (c) and Figs. [19] (c), an asymmetry 
in the polymeric stress between the upstream and down- 
stream regions is apparent. As seen in Fig. [T5](c), the 
asymmetry in the polymeric stress between the upstream 
and downstream regions develops gradually, while such 
asymmetry developments do not appear in the velocity 
and strain-rate field, as shown in Figs. [15] (a) and (b). 
The Weissenberg number We is larger than unity in the 
vicinity of the circle as shown in Fig. [18] (b) by thick 
lines. A nonlinear relationship is observed between D xy 
and a xy in the regions where We > 1 and even when 
We < 1. 

To highlight the nonlinear relationship between the 
strain rate D xy and the shear stress a xy in the asym- 
metric stress field shown in Fig. [15] (c) , we consider the 
deviation of the polymeric shear stress a xy from the New- 
tonian stress estimated using the zero-shear-rate viscosity 
of the polymer melt cr° = rj°D xy , as follows: 



Aa xy 



■xy\ 



(46) 



Figure [5D] shows the cases of (a) Aa xy > and (b) 
Aa xy < at t = 800 [tj. When the stress deviation 
Aa xy (r) is positive at a certain position r, the stress ex- 
hibits an overshoot coming from the elastic contribution 
of unoriented but well-extended polymer chains. When 
Aa xy (r) is negative, on the other hand, the stress repre- 
sents a shear-thinning behaviour. As seen in Fig. [2H1 (a) . 
stress overshoots are clearly observed in the upstream re- 
gions; these are indicated by the white dashed circles in 
the figure. In the regions where We > 1, the orientations 
of polymer chains cannot recover the isotropy. In the re- 
gions where We > 1, the inside of the solid lines in Fig. 
I2TH (b), we can clearly see the shear-thinning behaviour 
arising from the nonlinear relationship between the shear- 
rate and shear stress, e.g., a xy cx D™ y (0 < n < 1). 



Shear-thinning behaviour is also observed in the regions 
where We < 1, indicated by the black dashed circles in 
Fig. [2Q] (b). In the regions where We < 1, the orien- 
tations of polymer chains are expected to be isotropic, 
and the stress tensor should thus be proportional to the 
strain-rate tensor. However, the stress tensor is not pro- 
portional to the strain-rate tensor even in the regions 
where We < 1. Instead, nonlinear behaviours, i.e., stress 
overshoot and shear-thinning behaviour, appear in the 
regions where We < 1 because of the memory effect re- 
lating to the extensional relaxation. 

To quantitatively clarify the asymmetry of |v|, \D xy \ 
and \<r xy \ with respect to the j/-axis shown in Fig. I17[ we 
define a measure of asymmetry in a field / as 



J2 x =-H/aT, y i-H/a 1 1 /(*.!/) I ~ 1/0,2/)! 
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where x and y are the lattice coordinates in —H/a < x < 
H/a and —H/a < y < H/a, respectively, x = —x. The 
maximum value of Asym(/) at the steady-state /^ Ax is 
defined as /^ Ax = f t>u \fMAx(t)\dt/ J t>u dt where t s is 
the time when the system reaches the steady-state and 
/maxW is the maximum value of / in the system at the 
time t. 

Figure |2"T1 shows the time evolutions of (a) Asym(|v|), 
(b) Asym(\D xy \), and (c) Asymdcr^l) represented with 
dotted line, broken line, and solid line, respectively, for 
the polymer melt (P) with Z — 7 and rjd/rf 1 = 0.1 and 
the Newtonian fluid (N) with rf = r/^. Transient be- 
haviour can be seen in Fig. [21] before t = 400[t o ] in the 
polymer melt and t = 200 [t c ] in the Newtonian fluid. Af- 
ter t — 400 [t c ] in the polymer-melt case (or t — 200 [t e ] 
in the Newtonian case), these asymmetry measures ran- 
domly oscillate around each mean value in the steady- 
state. Thus, the time in steady-state t s is estimated to 
be 400 [t c ] in the polymer- melt case and 200 [t e ] in the 
Newtonian case. The asymmetry measure of \<J xy \ of the 
polymer melt shows a typical time-evolution behaviour, 
i.e., the flow history develops an asymmetry in the stress 
field. As shown in Fig. 1211 the asymmetry measure of 
\a xy \ of the polymer melt is much higher than that of 
|v| and \D xy \ in the steady-state. Reflecting the larger 
noise in Fig. [15] (b) than in Fig. [T5] (a), the asymmetry 
measure of \D xy \ of the polymer melt is slightly higher 
than that of |v| in the steady-state. The noise of the data 
raises the asymmetry measure. Even after accounting for 
the noise of the data, the asymmetry measure of \cr xy \ of 
the polymer melt is much higher than that of the others. 
The asymmetry measures of the Newtonian fluid at the 
steady-state have small values as shown in Fig. [31] (N) . 
This small asymmetry is caused by the asymmetric dis- 
tribution of the positions of fluid particles constituting 
the Lagrangian fluid. The asymmetry measures of |v| 
and \D xy \ of the polymer melt are slightly higher than 
those of the Newtonian fluid. The velocity and strain 
rate do not directly reflect memory of the flow history, 
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FIG. 18: Transient flow behaviour in the viscoelastic polymer melt. Contour maps of magnitudes of (a) velocity |v|[a/i e ], (b) 
strain rate lA^Kl/ie], and (c) shear stress over zero- shear viscosity \o xy \/rf[l/t c ] at 100, 200, 400, and 800 [t c ] are presented. 
Contour levels of these figures are 0.02 (solid line), 0.03 (broken line), 0.05 (dotted line), and 0.07 (dash-dot line) [a/t c ] in (a), 
and 0.002 (solid line), 0.004 (broken line), 0.006 (dotted line) [l/t e ] in (b) and (c). The thick solid lines in (b) indicate We = 1. 



although the memory indirectly affects them through the 
macroscopic stress, which is influenced by the memory of 
the strain rate. Therefore, the asymmetry measures of 
both |v| and \D xy \ are weaker than that of \u X y\ in the 
polymer-melt case. 



The asymmetry in the stress distribution is caused by 
the strain-rate- history-dependent stress, i.e., the mem- 
ory effect^ In Ref. [l?], we discussed the relationship 
between the stress relaxation time of a polymer melt in 
the homogeneous bulk and the retardation time of the 
stress field from the strain-rate field. We found that the 
retardation time is comparable to the stress relaxation 
time that actually corresponds to the Rouse relaxation 
time, i.e., the extensional relaxation time in the entan- 
gled polymer melt. 



V. SUMMARY AND FUTURE PERSPECTIVES 

We have developed multiscale bridging simulations for 
simple liquids, one-dimensional polymeric liquids, and 
two-dimensional polymeric liquids. In our multiscale 
modelling, the macroscopic flow behaviours are calcu- 
lated using the CFD method. However, instead of using 
any constitutive equations, a local stress is generated us- 
ing a microscopic (or mesoscopic) simulation, e.g., MD 
simulation and coarse-grained polymer dynamics, asso- 
ciated with a local fluid element according to local flow 
variables obtained by the CFD calculation. This type of 
multiscale modelling is expected to be useful in predict- 
ing the macroscopic behaviours of complex fluids such as 
polymeric liquids and other complex softmatters in which 
the constitutive relations are unknown. 

A multiscale simulation method of MD and CFD for 
simple fluids is developed in Section II. In this method, 
a usual finite- volume scheme with a lattice mesh is used 
for the CFD level, and each lattice is associated with a 
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FIG. 19: Transient flow behaviour in the Newtonian fluid. Contour maps of magnitudes of (a) velocity |v|[a/t e ], (b) strain rate 
|£) B j,|[l/te], and (c) shear stress over zero-shear viscosity Icr^j, l/r? [l/t e ] at 100, 200, 400, and 800 [t e ] are presented. In case of 
the Newtonian fluid, Figures (b) and (c) are abbreviated to the same figures because \D xy \ = \<J xy \/rf . Contour levels of these 
figures are 0.02 (solid line), 0.03 (broken line), 0.05 (dotted line) [a/t e ] in (a), and 0.002 (solid line), 0.004 (broken line), 0.006 
(dotted line) [l/i e ] in (b) and (c). 




FIG. 20: Greyscale contour maps of Aa xy = \a xy \ — \u l l y \[a c \ at t — 800[t c ] are separately presented for cases in which the values 
of Act xy are (a) positive and (b) negative. The regions surrounded by the solid lines (white in (a) and black in (b)) represent 
the regions in which We > 1. The dashed circles (white in (a) and black in (b)) indicate regions where the stress overshoot 
and shear-thinning behaviour caused by the memory effect in the stress are clearly observed in (a) and (b), respectively. 



small MD cell that generates a "local stress" according 
to a "local flow field" obtained from the CFD calcula- 
tion. The bridging scheme of CFD and MD is based on 
the local equilibrium assumption without memory (i.e., 
a local stress immediately attains the steady state after 
a short transient time according to a given flow field) 
[Fig. (3J. The cavity flows of a Lennard- Jones liquid in 
a square box are calculated for large efficiency factors, 
e.g., Ax//MD=At/TMD=8, and the results are compared 



with those of the Newtonian fluid and fluctuating hy- 
drodynamics. With these efficiency factors, multiscale 
simulations are performed 8 2 x 8 times faster than the 
full MD simulation. It is demonstrated that, although 
the instantaneous profiles of macroscopic quantities in- 
volve large fluctuations for large efficiency factors, the 
fluctuations are perfectly smoothed out by taking the 
time averages [Figs. [6] and [7]. The spectrum analyses of 
fluctuations are carried out, and the fluctuations arising 
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FIG. 21: Time evolution of asymmetry measure of (a) |v| (dotted line), (b) \D xy \ (broken line), and (c) \<J xy \ (solid line) for 
(P) polymer melt flow with Z = 7 and r/d/?? = 0.1, and (N) Newtonian flow with r; = rjd- 



in the multiscale simulations are found to be consistent 
with those of the fluctuating hydrodynamics [Figs. [8] and 
[9] . We also found a relation between the fluctuation and 
efficiency factor as Eq. (|22[) . 

In Sec. Ill, the multiscale simulation method for sim- 
ple fluids is extended to one-dimensional flows of poly- 
meric liquids in parallel plates, in which the macroscopic 
quantities are assumed to be uniform in the flow direc- 
tion parallel to the plates. This assumption allows us to 
not trace the advection of memory of polymer configu- 
rations in each fluid element but instead to calculate the 
coupling of macroscopic flow behaviours to the micro- 
scopic dynamics of polymer chains using the usual fixed- 
mesh system at the CFD level. As in the simple fluids, 
a local stress is calculated by the non-equilibrium MD 
simulation in a small MD cell associated with a mesh 
interval of the CFD calculation. However, the MD sim- 
ulations are performed in the time-step duration of the 
CFD calculation at each time step of CFD, and the fi- 
nal configuration of polymers obtained in each MD cell 
is retained as the initial configuration in each MD cell at 
the next time step. Thus, the memory of configurations 
of polymer chains is traced at the MD level so that the 
memory effects are correctly reproduced [Fig. 110] . The 
behaviour of a glassy polymer melt in rapidly oscillating 
plates is calculated for various amplitudes and oscillation 
frequencies, and the velocity profiles and local Theologi- 
cal properties of the melt are investigated. The results 
of the multiscale simulations for a plate with an oscil- 
lation of small amplitude are also compared with those 
of the analytical solution for a plate with an infinitesi- 
mally small amplitude oscillation so as to demonstrate 
the validity of our multiscale simulations. The velocity 
profiles become thinner as the amplitude of oscillation of 
the plates increases due to the shear thinning occurring 
near the oscillating plate, where the local strain rates are 
large [Figs. QT] and [12]. It is also demonstrated that 
the velocity profiles of the melt obtained by the multi- 
scale simulations approach that of the linear analytical 
solution as the amplitude of oscillation of the plates de- 



creases [Figs. [T3l and [14] . The local rheological prop- 
erties are investigated in terms of the complex viscosity 
[Fig. [15] . The local complex viscosities are demonstrated 
to become uniform and to approach those of the linear 
response as the amplitude of oscillation decreases. How- 
ever, for large amplitudes, the local rheological properties 
are quite spatially varied. It is found that at high oscilla- 
tion frequency and large amplitudes of oscillation of the 
plates, the melt has three different rheological regimes, 
i.e., viscous liquid, viscoelastic liquid, and viscoelastic 
solid regimes. 

In Sec. IV, the multiscale simulation method, com- 
posed of a Lagrangian fluid dynamics simulation and a 
coarse-grained polymer simulation, is presented and ap- 
plied to the entangled polymer-melt flow, in which the 
advection of microscopic variables affects its flow be- 
haviour. The Lagrangian fluid dynamics simulation en- 
ables us to treat the strain-rate-history-dependent stress 
field that is evaluated by a large number of microscopic 
simulations. A local stress of the entangled polymer melt 
is obtained by a slip-link model that can efficiently calcu- 
late the dynamic response of the entangled polymer melt 
to an imposed strain rate in the homogeneous bulk. As a 
demonstration, we consider the transient flow of the poly- 
mer melt passing a circular object in a two-dimensional 
periodic system. In order to use the Lagrangian fluid 
dynamics simulation at the macroscopic level, we can 
trace the polymer conformations in the inhomogeneous 
flow field, and show that the effect of memory that is 
preserved as the polymer conformations in a fluid par- 
ticle appears as a macroscopic stress distribution with 
an asymmetry between the upstream and downstream 
regions around the circular object [Fig. 115] . Nonlinear 
behaviours, represented by stress overshoot and shear- 
thinning behaviours, are observed in the regions where 
We < 1 because of the memory effect relating to the 
extensional relaxation [Fig. [20]. To quantify the asym- 
metry between the upstream and downstream regions 
around the circular object, we define an asymmetry mea- 
sure and investigate its transient behaviour. The asym- 
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metry measure of the polymeric stress is found to show 
typical transient behaviours: i) the asymmetry measure 
of the stress increases before t = t s and ii) the asymme- 
try measure approaches a steady-state value after t = t s , 
while the asymmetry measure of the velocity and the 
strain rate do not represent transient behaviour [Fig. [21] . 

In the present work, we have developed a multiscale 
simulation of the multi-dimensional flows of polymeric 
liquids involving the advection of memory. Tracing of 
the memory and its advection is quite important not only 
in polymeric flows but also in many kinds of softmatter 
flows. In this work, we discuss only isothermal fluids with 
a uniform density. We have also assumed the non-slip 
boundary condition for the macroscopic velocity at the 
wall. Thus, there remain two important directions of de- 
velopment of multiscale modelling for the complex flows 
of softmatters. One is the treatment of the couplings of 
various macroscopic quantities, e.g., temperature, den- 



sity, and velocity. In particular, the inclusion of the 
temperature variation may be important for polymeric 
liquids. The other important developmental direction 
involves consideration of the role of microscopic dynam- 
ics at the interface in determining macroscopic flow be- 
haviours. Some simulation methods have been proposed 
for this purpose.— ~— Incorporating one of these into our 
multiscale simulation method may enable us to treat slip 
motion, adhesion, and anchoring at the interface together 
with the macroscopic flow behaviours. Furthermore, al- 
though the high computational expense of describing 
three-dimensional systems prevents us from simulating 
the three-dimensional flows, the presented multiscale for- 
mulations allow consideration of three-dimensional flows 
of soft materials in principle. In the near future, develop- 
ments in computer architecture and more parallelisation 
should enable us to account for such dense multiscale 
systems. 
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